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.^ . Abstract Coupled flux transport and magncto-frictional simulations are ex- 

tended to simulate the continuous magnetic field evolution in the global solar 
corona for over 15 years, from the start of Solar Cycle 23 in 1996. By simplifying 
the dynamics, our model follows the build-up and transport of electric currents 
P^ ■ and free magnetic energy in the corona, offering an insight into the magnetic 

structure and topology that extrapolation-based models can not. To enable 
(-h . these extended simulations, we have implemented a more efficient numerical 

fS || grid, and have carefully calibrated the surface flux transport model to reproduce 

the observed large-scale photosphcric radial magnetic field, using emerging active 
j^ ' regions determined from observed line-of-sight magnctograms. This calibration 

^ ■ is described in some detail. In agreement with previous authors, we find that the 

standard flux transport model is insufficient to simultaneously reproduce the 
observed polar fields and butterfly diagram during Cycle 23, and that additional 
effects must be added. For the best-fit model, we use automated techniques to 
detect the latitude-time profile of flux ropes and their ejections over the full 
solar cycle. Overall, flux ropes are more prevalent outside of active latitudes 
but those at active latitudes are more frequently ejected. Future possibilities for 
f*^ | space-weather prediction with this approach are briefly assessed. 

Keywords: Coronal mass ejections, theory; Magnetic fields, corona; Magnetic 
fields, models; Magnetic fields, photosphere; Solar cycle, models 



1. Introduction 



?—i ■ Modelling the Sun's coronal magnetic field over the full solar cycle is important 

because it acts as a driver of space-weather events, and changes significantly 
over the cycle, as well as from one cycle to the next. It is fundamentally time- 
dependent, as manifested in the variations of almost all observed properties of 
the Sun, including sunspot number, magnetic flux, rates of flares and coronal 
mass ejections (CMEs), and even the total solar irradiance (Willson and Hudson, 
1991). The coronal magnetic field is driven by activity in the solar interior, 
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including large-scale flows and convection in the photosphere as well as the 
periodic emergence of new active regions. 

Previous studies of the coronal magnetic evolution over months to years have 
mostly used potential-field extrapolations (Altschulcr and Ncwkirk, 1969; Schat- 
ten, Wilcox, and Ness, 1969). From the space-weather viewpoint, these models 
give a first approximation of the Sun's open flux (Wang and Sheeley, 2002) 
and the heliospheric current sheet (Hoeksema, Wilcox, and Scherrer, 1983), and 
have been coupled to models of the heliosphere (Arge and Pizzo, 2000; Luhmann 
et at, 2002; Pizzo et at, 2011). They have also been applied to predict CME 
rates using the topology of magnetic nulls and "breakout" configurations (Cook, 
Mackay, and Nandy, 2009). However, their value for predicting flares and CMEs 
is limited since they allow no free energy, and the lack of electric currents leads 
them to underestimate the open flux, particularly during active periods (Riley, 
2007; Yeates et at, 2010a). 

Attempts to move beyond potential-field models and extrapolate more gen- 
eral magnetic equilibria suffer from the problem of non-uniqueness of magnetic 
topology for a given distribution of observed magnetic field in the photosphere. 
For example, the current-sheet source-surface (CSSS) model can better match 
the observed open flux by allowing certain forms of electric currents (Zhao and 
Hoeksema, 1995; Jiang et at, 2010a), but these are chosen for mathematical 
convenience rather than any particular physical basis. On the other hand, when 
trying to extrapolate nonlinear force- free fields from photospheric vector magne- 
tograms, different numerical extrapolation codes tend to produce different results 
(De Rosa et at, 2009). These difficulties have led to a more pragmatic approach 
for nonlinear force-free modelling of local structures such as coronal cavities or 
sigmoids, whereby the computation is initialised with a flux rope structure in 
the corona (van Ballcgooijcn, 2004; Su et at, 2009; Savcheva, van Ballcgooijen, 
and DcLuca, 2012). 

In recent years, global magnetohydrodynamic (MHD) models that more re- 
alistically account for thermodynamic properties of the plasma have become 
practical, allowing for comparison with observed emission at various wavelengths 
(Lionello, Linker, and Mikic, 2009; Rusin et at, 2010; Downs et at, 2010). 
However, it is not yet practical to simulate the temporal evolution of the corona 
for many months, and MHD models are generally limited to finding individual 
equilibria by relaxing from an initial condition appropriate for a given day. So 
far, the initial conditions for the magnetic field have been potential-field extrap- 
olations. Thus the MHD models inherit the topology of the potential field, and 
do not account for the gradual build up of the magnetic topology over time, or 
any long-term memory of previous interactions that remains imprinted in the 
coronal field. 

The approach of our model is to gain new insight into the magnetic structure 
and topology by not just extrapolating from photospheric data at a single time, 
but simulating the coronal field in a time-dependent way, using a simplified 
approximation to the real coronal evolution. The technique was originally in- 
troduced to study the formation of filaments (van Ballcgooijcn, Priest, and 
Mackay, 2000; Mackay, Gaizauskas, and van Ballcgooijcn, 2000; Mackay and 
van Ballcgooijcn, 2001, 2005, 2006), whose magnetic structure is non-potential 
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and conjectured to depend on the build-up and transport of coronal magnetic 
hclicity over time (van Ballegooijen and Martens, 1989). It was extended to 
model the global corona by Yeates, Mackay, and van Ballegooijen (2008). The 
model effectively produces a continuous sequence of nonlinear force- free fields, in 
response to flux emergence and shearing by photospheric footpoint motions. As 
such, a particular magnetic topology is chosen automatically at each time step 
by the evolutionary history, thus providing a physically motivated solution to 
the non- uniqueness problem. Moreover, the model - henceforth the NP ("non- 
potential") model - offers interesting possibilities for modelling and predicting 
space weather, because it allows for the build up and transport of free magnetic 
energy, electric currents, and magnetic helicity. In the model - as is hypothesized 
in the real corona - helicity tends to concentrate in flux rope structures overlying 
photospheric polarity-inversion lines. When too much hclicity accumulates, the 
flux ropes "erupt" and are ejected through the outer boundary of the simulation 
domain (Mackay and van Ballegooijen, 2006; Yeates and Mackay, 2009). 

The developments in this article are threefold. Firstly, we extend the sim- 
ulations to a continuous 15 year evolution (Section 2). Our previous study of 
solar cycle variations in the NP model was limited to six separate six-month 
runs (Yeates et al., 2010b). The new simulation allows for longer-term magnetic 
memory to take effect; our initial finding that this affects the chirality of high- 
latitude filaments has been described elsewhere (Yeates and Mackay, 2012). 
Secondly, we describe how the surface flux transport component of the model 
- the lower boundary condition - must be carefully calibrated to observations 
for such long simulations (Section 3). Thirdly, we consider the distribution of 
flux ropes and flux rope ejections over the 15-year simulation (Section 4). We 
conclude in Section 5 with a brief discussion of the future prospects for space 
weather forecasting using this approach. 



2. Coronal Magnetic Model 

The non-potential (NP) model couples surface flux transport to magncto-frictional| 
relaxation in the overlying corona (van Ballegooijen, Priest, and Mackay, 2000). 

2.1. Formulation 

The large-scale mean coronal magnetic field [Bo = V x Ao] is evolved by the 
induction equation 

^=v xB„-E , (1) 

where we neglect ohmic diffusion and the mean electromotive force [En] describes 
the effect of unresolved small-scale fluctuations. Following van Ballegooijen and 
Cranmcr (2008), we apply a hypcrdiffusion 

E o = -S V -(^BoVao), (2) 
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where 

Bo -jo / n 

a a = 2 (3) 

is the current helicity density, with jo = V x Bo the current density, and 7/4 = 
10 11 km s _1 . This form of hyperdiffusion preserves magnetic-helicity density 
[A • B ] in the volume and describes the tendency of the magnetic field to 
relax to a state of constant «o (Boozer, 1986; Bhattacharjee and Hameiri, 1986), 
although such a state is never reached in the global simulation. The velocity is 
determined by the magneto-frictional technique (Yang, Sturrock, and Antiochos, 
1986; Craig and Sneyd, 1986) as 



1 jo x_B 
-•o 



v = ^5 !" v out(r)e r . (4) 



This replaces the full momentum equation and allows numerical solution of the 
model over months and years. The first term enforces relaxation towards a force- 
free equilibrium [jo x Bo = 0]. The second term is a radial outflow imposed only 
near the outer boundary [r = 2.5R0] to represent the effect of the solar wind 
radially distending magnetic field lines (Mackay and van Ballegooijen, 2006). 

2.2. Photospheric Boundary Condition 

On the photospheric boundary r = Rq, the magneto-frictional velocity is not 
applied, and instead the radial magnetic field [Bq t ] is evolved by the surface flux 
transport model (Sheeley, 2005; Mackay and Yeates, 2012). In spherical polar 
coordinates (r, 0, 4>), the vector potential evolves according to 

-—— = u^Bor - +Se(9,4>,t), (5 

at R© sm v ocp 

— = - Ue B Qr + —^ F + S 4> {e,cp,t). (6) 

Here D is a (constant) diffusivity modelling the random walk of magnetic flux 
owing to the changing supcrgranular convection pattern (Lcighton, 1964). In 
Section 3, we experiment with an additional exponential decay of B$ r (Schrijver, 
De Rosa, and Title, 2002). The differential rotation velocity u^ = 17(0)R Q sin# 
uses the observationally-determined Snodgrass (1983) profile 

0,(6) = 0.18 -2.3 cos 2 6- 1.62 cos 4 6>deg day" 1 , (7) 

written in the Carrington frame. For the basic meridional flow we assume the 
form of Schiisslcr and Baumann (2006), namely 

16 

MO) = «o yy75 sin ( 2A ) cx p ( n - 2 I A D . ( 8 ) 

where A = n/2 — 9 is latitude and uq is a constant controlling the flow amplitude. 
In Section 3, we experiment with activity-dependent perturbations of this flow 
profile. 
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The source terms Sg and S^ in Equations (5) and (6) represent the emergence 
of new active regions, and are necessary to maintain an accurate description of 
the observed surface Bq t over the continuous 15 year simulation. Rather than 
specify a functional form for these source terms, we insert individual bipolar 
magnetic regions with properties chosen to match those in observed synoptic 
magnctograms (Yeates, Mackay, and van Ballcgooijen, 2007). The inserted bipo- 
lar regions take an idealised three-dimensional form (Yeates, Mackay, and van 
Ballegooijen, 2008). In the future, we hope to incorporate more detailed models of 
the structure of individual active regions, built up in a time-dependent manner. 
For the simulations described here, we use synoptic normal-component magne- 
tograms from US National Solar Observatory, Kitt Peak. Until 2003, these were 
taken with the older Vacuum Telescope, and from 2003 onward with Synoptic 
Optical Long-term Investigations of the Sun (SOLIS). We insert a total of 2040 
bipoles between Carrington Rotation CR1911 (June 1996) and CR2122 (April 
2012). We do not insert any bipoles to replace those missed during the three 
data gaps (CR2015-16, CR2040-41, CR2091). 

The 2040 bipolar regions are summarised in Figure 1. Figure 1(a) shows their 
latitude-time distribution, along with their leading/following polarity. The ma- 
jority polarity reverses with each 11-year cycle, as it should. Figure 1(b) shows 
the range of sizes and fluxes of bipoles in our dataset. Our semi-automated tech- 
nique identifies bipoles using only flux above 50 gauss, thus imposing a lower size 
cut-off. Figure 1(c) shows the latitude distribution of bipole tilt angles [7] (the 
angle with respect to the Equator of the line connecting leading and following 
polarity centroids). The solid line shows a linear fit 7 = —0.41° + 0.459A, which 
is consistent with magnetogram observations of Wang and Sheeley (1989), who 
used NSO/KP data for Cycle 21, or Stenflo and Kosovichcv (2012), who used 
SOHO/MDI magnctograms. White-light studies of tilt angles return a rather 
lower slope (Dasi-Espuig et al., 2010). We feel that the distribution derived 
from magnctograms should be more appropriate for our application, although 
we find in Section 3 that reducing the tilt angles by 20 % is a straightforward way 
to calibrate the flux-transport simulation without the need for more complex 
effects (Cameron et al., 2010). For the 20% reduced tilt angles, the linear fit is 
7 = -0.33°+0.367A. 

2.3. Numerical Methods 

The equations are solved for the vector potential [Ao] on a staggered grid, using 
a flux-conserving (finite-volume) scheme for the advection terms. To avoid the 
problem of grid convergence near the poles, we have newly incorporated a vari- 
able grid (see Appendix A for details). This greatly reduces the computational 
time. For the simulations described in this article, we use a resolution of 192 cells 
in longitude at the Equator, corresponding to an angular resolution of 1.875°, 
and 28 cells in radius. At the polar grid boundary (6 ps 0.67°), there are only 12 
cells in longitude, owing to the variable grid. Boundary conditions are described 
in Appendix B. The coronal magnetic field is initialised using a potential-field 
extrapolation for 15 June 1996, derived from the synoptic magnetogram for 
CR1910 (Yeates, Mackay, and van Ballegooijen, 2007). The code is parallelised 
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Figure 1. Properties of the 2040 magnetic bipolcs determined from NSO/KP synoptic mag- 
nctograms: (a) locations in time and latitude, with colour/symbol showing polarity and symbol 
size proportional to flux; (b) bipole flux against bipole size; (c) tilt angles against latitude. In 
(c), the solid line is a linear fit to the measured tilt angles. 



with OpenMPI and the full 5822-day run took approximately five days with 48 
cores. 



3. Observational Calibration 

Since we do not reset the photospheric magnetic field to observed magnetograms 
once the simulation has begun, it is important to calibrate the surface flux- 
transport model to reproduce the observed long-term evolution on the Sun. 
This is particularly true for simulations as long as 15 years, so we look at this 
issue in detail here. Unfortunately, it is a delicate balance between the proper- 
ties of newly-emerging active regions and the chosen profiles for the transport 
processes of meridional flow and supergranular diffusion. Differential rotation is 
better constrained by observations and we do not consider its variation. Pre- 
vious parameter studies of the flux-transport model have been undertaken by 
Baumann et al. (2004), and for Cycle 23 by Schrijver and Liu (2008) and Jiang 
et at (2011), although here we are additionally constrained by the individual 
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measured properties of the bipolar regions. Since there is no feedback of the coro- 
nal (magneto-frictional) component of the model on the surface flux-transport 
component, the calibration can be done by running only the latter. 

To calibrate the surface simulations we use several numerical indicators for 
the photosphcric magnetic field: 

(i) Cross-correlation of total (unsigned) magnetic flux [<j> ] with the observed 
time series. 

(ii) Least value of $0 around the Cycle 23 Minimum (ca. 2009). 

(iii) Cross-correlation of axial dipole strength B ax = 6 10 with the observed time 
series. 

(iv) Time of sign-reversal in the axial dipole strength. 

(v) Peak negative value of the axial dipole strength during Cycle 23. 



(vi) Cross-correlation of equatorial dipole strength B cq = -\/|&i,i| 2 + | t>i, 1 1 2 = 

a/^I&i.iI with the observed time series. 

(vii) Cross-correlation of longitude-averaged field (Bo r ) = j Q BorR© sin 6 dcf> 
with the observed butterfly diagram. 

Here bi, m are the complex- valued spherical-harmonic components of Bo r . Be- 
cause of the smoothing effect of supergranular diffusion, the total simulated 
flux [Bor] tends to be lower than the total magnetogram flux. To account for 
this, our cross-correlations (a) and (g) use a smoothed version of the observed 
magnetogram. obtained by removing spherical harmonic components above order 
& = f 28. Higher-resolution simulations in future will need to incorporate a more 
realistic random walk model for discrete flux concentrations (e.g., Wordcn and 
Harvey, 2000; Schrijver, 200f). To facilitate cross-correlation of time series in 
(a), (c), and (f), the simulation data are first interpolated to the same times 
as the magnetogram observations. For cross-correlation of the two-dimensional 
butterfly diagrams, the simulation data are first interpolated to the same times 
and latitudes as the observations. 

Table 1 lists values of the numerical indicators for a representative set of 
test runs, which are described in more detail below. For comparison, the total 
unsigned flux [$©], axial dipole strength [-B a x], and equatorial dipole strength 
[B eq \ for each run are shown in Figure 2. Butterfly diagrams of (B$ r ) in each 
run are shown in Figure 3. 

3.f . Reference Case 

Consider first the reference case Run 0, with peak meridional flow uq = f f ms _1 
(consistent with Cycle 23 observations by Hathaway and Rightmire, 20f 0) and 
diffusivity D = 450 km s _1 . The time series for Run are shown by black 
curves in all panels of Figure 2, while corresponding values for the observed 
magnetograms are shown by thick grey lines. It is clear that $© is too high from 
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Table 1. Performance of surface flux 


transport calibration 


runs. 




Run 

(i) (ii) [10 23 Mx] 


Indicators 
(iii) (iv) [years] 


(v) [gauss] 


(vi) (vii) 


Observed 1.520 
Smoothed 0.895 


2000.04 


-5.47 


0.945 






0.941 


1.314 


0.934 


1999.71 


-14.49 


0.688 


0.693 


A35 


0.965 


0.392 


0.938 


2001.19 


-6.47 


0.570 


0.590 


B200 


0.968 


1.146 


0.951 


1999.93 


-9.85 


0.684 


0.709 


C80 


0.960 


1.119 


0.933 


2000.08 


-10.32 


0.709 


0.711 


D80 


0.954 


1.147 


0.931 


2000.08 


-10.88 


0.687 


0.708 


E 


0.967 


0.898 


0.898 


2001.04 


-8.36 


0.673 


0.695 


F5 


0.966 


0.553 


0.913 


1999.19 


-9.74 


0.682 


0.706 


C80F10 


0.967 


0.702 


0.966 


1999.64 


-8.65 


0.706 


0.726 



about 2004 onwards, during the declining phase. This is clearly caused by an 
excess in B ax rather than B cq . Indeed, the butterfly diagram (Figure 3) shows 
that Run builds up unrealistically large polar caps. The low polar field of Cycle 
23 is not reproduced by this standard model. A second inconsistency in Run is a 
deficit in the equatorial dipole strength B cq around 2002-2004. This is present in 
all of the calibration runs and appears to be due to an underestimate of emerged 
flux in a handful of particular active regions near the Equator, to which B oq is 
sensitive. This is likely caused by our simplified method of extracting bipolar 
regions, but fortunately does not have a significant effect on $ Q or on the later 
polar field produced. 

3.2. Meridional Flow Speed 

A possible solution to the excess polar field problem is to retain the standard flux 
transport model but to alter the meridional flow speed [uq] or diffusivity [D] . Dif- 
ferential rotation is better constrained and does not significantly affect the dipole 
strengths (Baumann et al., 2004). In order to reduce the polar field production, 
one must increase uo, so as to reduce cancellation at the Equator (DcVorc, 
Shcclcy, and Boris, 1984). In Run A35, the required speed of u = 35ms -1 is 
used to bring B ax down to the observed level (Figure 2). But in addition to being 
inconsistent with observations (Hathaway and Rightmire, 2010), this causes a 
very late reversal of B ax , and too low a level of B cq over much of the cycle. The 
resulting butterfly diagram is also poorly correlated with observations (Figure 
3). 

Schrijver and Liu (2008) and Jiang et al. (2011) were able to improve the 
match of their simulations to observed dipole moments/polar fields by changing 
the meridional flow. However, when we run a simulation with our measured 
bipolc properties and the diffusivity (250 km s -1 ) and meridional flow profile of 
Jiang et al. (2011), we find that the axial dipole strength is close to that of Run 0, 
even with their 55 % enhanced flow speed (17ms -1 ). Moreover, correlation with 
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Figure 2. Time series of total unsigned photospheric flux [$0], axial dipolc component 
[Sax], and equatorial dipolc component [B eq ] for the various surface calibration runs. Ob- 
served (NSO/KP) values are shown by thick grey lines (in the top row, the dashed grey 
line shows value from the original magnetograms and the solid grey line that from smoothed 
magnetograms) . 



the observed butterfly diagram is again poor. If we implement the meridional- 
flow profile of Schrijver and Liu (2008), with a stronger latitudinal gradient in 
ug at the Equator, and uq = 15ms -1 , we still obtain similar results. Note that 
neither of those studies used observations of individual bipolar magnetic regions, 
as we do here. 

3.3. Supcrgranular Diffusivity 

In Run B200, we instead reduce D to 200 km s _1 . This produces slightly high 
$0 compared with the observations throughout the cycle, although the butterfly 
diagram and B cq are a better fit to the observations. In the declining phase, S ax 
is still rather high. Here we note a discrepancy with the parameter study by 
Baumann et al. (2004) for our u$ profile and bipolc properties: as we increase 
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Figure 3. Magnetic butterfly diagrams showing the longitude-average (Bq t ) for each calibra- 
tion run, and for the observed KP magnctograms (bottom right). In each case the greyscalc 
saturates at ±6 G. Data gaps and problems with high-latitude measurements are evident in 
the observed butterfly diagram. 
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D, the peak axial dipole B ax increases, whereas they found the maximum polar 
field strength to decrease over the same range of D. 

3.4. Bipole Properties 

Although our bipolar regions are constrained individually by observations, we 
also consider systematically modifying their properties. One can reduce the polar 
field either by reducing the bipole fluxes, or by reducing the average tilt angle 
so that each region contributes less to the axial dipole (Wang, Sheeley, and 
Lean, 2000). In Run C80, all tilt angles were reduced by 20%, while in Run 
D80 all bipole fluxes were reduced by 20%. Both runs produce a comparable 
improvement in B ax , with accurate reversal times. Overall, Run C80 produces 
better results than Run D80, because the latter causes too much reduction in B eq 
and consequently in $ Q , particularly during active periods. Jiang et al. (2011) 
found that a 28 % decrease in their tilt angles (as compared to previous cycles) 
produced a reasonable polar field in Cycle 23. However, even in Run C80 the 
peak of -B ax remains too strong compared to the observations, so we are led to 
consider an additional change to the model. 

3.5. Exponential Flux Decay 

In run F5, we try adding additional exponential decay terms of the form — Aqq/t 
and — Ao^/t to Equations (5) and (6) respectively, leading to an exponential 
decay on a timcscale r — 5 years. Such a decay has previously been introduced 
by Schrijver, Dc Rosa, and Title (2002) in flux-transport simulations of the 
past 340 years, in order to maintain regular polar-field reversals when cycles 
vary in strength from one to the next. Baumann, Schmitt, and Schiisslcr (2006) 
have introduced a similar enhanced decay, proposing the physical explanation to 
be volume diffusion of the surface magnetic field in the three-dimensional solar 
interior. Figure 2 shows that this decay is able to bring B ax down to the observed 
level in 2010 while maintaining a good correlation with the observed butterfly 
diagram, although it has the side-effect of bringing the reversal time of i? ax too 
early. By combining a weaker exponential decay term (r = 10 years) with 20 % 
reduced tilt angles (as in Run C80), Run C80F10 produces a better compromise 
and an even stronger correlation with the observed butterfly diagram. This is 
the run chosen for driving the coronal simulations in Section 4, and used by 
Yeates and Mackay (2012). 

3.6. Time- varying Meridional Flow 

An alternative improvement to the flux-transport model is to introduce tem- 
poral and spatial fluctuations in the meridional flow. Either a high-latitude 
countercell (Jiang et al., 2009) or converging flows towards active regions (Jiang 
et at, 2010b) can reduce the resulting polar field. Although we do not carry 
out an exhaustive investigation here, we illustrate (Run E) the effect of such 
converging flows using the method of Cameron and Schiisslcr (2010). An ax- 
isymmctric perturbation is added to the steady flow uq to model the net effect 
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of non-axisymmetric perturbations (De Rosa and Schrijver, 2006), and takes the 
form 

u'(\,t) = co(±(\B 0r \)(\,t)\, (9) 

where (|i?or|)(A,i) is the longitudinal average of |Bor|; to which we also apply a 
Gaussian smoothing in latitude at each time. Following Cameron and Schiissler 
(2010) we set cq = 10 ms _1 gauss _1 deg. The Gaussian smoothing is implemented 
(at each time step of the main simulation) by taking 80 one-dimensional (in A) 
diffusive steps with diffusion coefficient 2.4 x 10 7 km 2 s _1 . Cameron and Schiissler 
(2010) show that this form of u'(X, t) leads naturally to variation of the Pi(cos6) 
component of the total meridional flow over the cycle. It thus reproduces both 
the observations of Hathaway and Rightmire (2010), who found a variation of 
the Py component using MDI tracking, and Basu and Antia (2010), who both 
substantiated the Hathaway and Rightmire result and found evidence for inflows 
towards active regions, using MDI helioscismology. For our bipole properties, 
Figure 2 shows that these flow perturbations lead to reasonable $q, and reason- 
able B ax later in the cycle. But the reversal of B ax is much later than observed, 
because poleward transport is reduced during the rising phase of the cycle. This 
leads also to poorer correlation with the observed butterfly diagram. Hence we 
do not pursue this route here, although temporal variations in meridional flow 
will be important to develop in the future, and particularly to take into account 
variations between different cycles (Wang, Lean, and Sheeley, 2002; Schrijver 
and Liu, 2008). 



4. Flux Ropes and Eruptions 

As an application of the NP model, we focus here on the formation and ejection 
of magnetic-flux ropes over the 15-year simulation, driven by the optimal flux 
transport Run C80F10. Twisted flux ropes form naturally in the simulation 
when surface motions concentrate magnetic helicity above polarity inversion lines 
in the photospheric field. Previously, we have developed automated techniques 
to detect them, and have analysed the effect of the various coronal simulation 
parameters on their formation and ejection (Yeates and Mackay, 2009). We have 
also undertaken a detailed comparison with CME source regions observed in 
the extreme ultraviolet (Yeates et al., 2010b). Yeates, Constable, and Martens 
(2010) looked at how the flux rope statistics differed in different phases of Cycle 
23, using six-month "snapshots". Here, we present the distribution of flux ropes 
and flux-rope ejections over the continuous 15-year simulation. This is intended 
to be an illustration of the model capabilities, rather than a detailed parameter 
study. In particular, the model is not yet suitable for prediction of individual 
space-weather events (see Section 5). 

4.1. Definition and Automated Identification 

Identifying magnetic-flux ropes in complex three-dimensional magnetic fields is 
a non-trivial task. We do this once per day during the course of the simulation, 
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Figure 4. Flux rope points are identified by an inward magnetic tension force [T r ] and an 
outward magnetic pressure force [Pa- 



using a slight modification of the automated technique described by Yeates and 
Mackay (2009). The technique works by first identifying flux-rope points on the 
numerical grid satisfying certain criteria, before grouping these points into flux 
ropes using a clustering algorithm. We have settled on the following criteria for 
defining flux-rope points. At each point on the computational grid, we compute 
the normalised vertical magnetic tension force and pressure gradient, 






T r = ^B -V ^ , P P = -^-(i3L). (10) 



A grid point [r% , 8j ,4>k) is then selected if it satisfies the following five conditions: 

Pr{n-i,6 h <l>k) < -0.4, (11) 

Prin+uOj^k) > 0.4, (12) 

Trin^Oj^k) > 0.4, (13) 

Trin+uOjrfk) < -0.4, (14) 

|jo-Bo| > a*Bl (15) 

where a* = 0.7 x 10~ 8 m -1 . Thus a flux rope is defined as a twisted structure 
where magnetic pressure acts outward from the flux rope axis and magnetic 
tension acts inward (Figure 4). As an illustration, Figure 5 shows four snapshots 
of the corona in Run C80F10, with red/orange magnetic-field lines traced from 
identified flux-rope points. 

To identify flux-rope ejections, we use the automated post-processing pro- 
cedure of Yeates and Mackay (2009). This selects those flux-rope points with 
vo r > 0.5 kms -1 in the magncto-frictional code and clusters them both spatially 
and temporally into separate ejection events. To be classed as separate events, 
clusters must be separated by at least five days in time and have at least eight 
points (on the grid resolution used here). This yields a list of times, locations, 
and sizes of flux rope ejections during the simulation. 

4.2. Latitude-Time Distribution 

At any one time, there are multiple flux ropes, of varying sizes, present in the 
simulation. To show how these vary over latitude and time, Figure 6 summarises 
the results of the automated detection routines applied to Run C80F10. 
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(a)25-Jul-1999 




(c) 06-Apr-2008 



Figure 5. Snapshots of the NP model at a sequence of times during Run C80F10, including (a) 
Cycle 23 Maximum prior to polar reversal, (b) the declining phase, (c) Cycle 23 Minimum, and 
(d) Cycle 24. In each case, blue field lines are traced down from the source surface r = 2.5Rq, 
while red field lines are traced from (a subset of) flux rope points. 



Figure 6(a) shows the distribution of flux ropes in time and latitude. The 
quantity plotted is the filling factor of flux ropes, namely the proportion of grid 
points at each latitude and time that are identified as flux rope points. There 
are two features to notice about the flux-rope distribution. Firstly, the latitude 
range of flux ropes reflects the extent of PILs on the solar surface, and varies 
over the cycle. The distribution of PILs reaches its broadest extent during the 
"rush-to-the-poles" in the run up to polar field reversal in 2000-2001, after 
which high-latitude flux ropes disappear again. The second feature to notice is 
that the filling factor is greater outside of active regions, either at high latitudes 
or during the Minimum period from 2007-2010 when there were few active 
regions present. This anti-correlation with the magnetic butterfly diagram fits 
the picture of flux ropes forming as a result of gradual transport and build-up 
of hclicity by surface motions. In fact, the number of individual flux ropes docs 
increase by a factor of about 1.5 during Solar Maximum (Yeates, Constable, and 
Martens, 2010), but this is because a greater number of smaller flux ropes are 
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Figure 6. Latitude— time distribution of (a) flux ropes and (b) flux— rope eruptions, and (c) 
time series of the flux— rope eruption rate. In (a), the colour scale shows flux rope filling 
factor, while in (b), the colour scale shows the fraction of flux rope points erupting in a given 
latitude-time bin. Bins with no flux rope points are coloured black. In (c), the observed CME 
rate from the CDAW catalogue divided by three is plotted in red (see text). 
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present in and around active regions. Larger ropes are found in decaying flux 
regions where helicity has had time to concentrate. 

This observed distribution of flux ropes shares some similarities with observed 
butterfly diagrams of solar filaments, which also overlie PILs (d'Azambuja and 
d'Azambuja, f948; Mouradian and Soru-Escaut, 1994). Namely, the overall lat- 
itudinal range and the increase in the number of filaments with solar activity. 
However, Mouradian and Soru-Escaut (1994) find that, around the Minimum of 
Cycle 21, filaments persisted only at higher latitudes and not near the Equator. 
This appears to be at odds with the simulated distribution of flux ropes in 
Figure 6(a) in the years 2008-2010, as seen in Figure 5(c). This requires further 
investigation. 

Figure 6(b) shows the proportion of flux-rope points found to be erupting, 
within each latitude-time bin. Empty bins are black. Notice that the pattern is 
predominantly the inverse of Figure 6(a), and now correlates with the magnetic 
butterfly diagram. So we find that, although less space is filled with flux ropes at 
active latitudes, the flux ropes that are present at active latitudes are more likely 
to erupt. Hence the eruption rate is modulated by solar activity (black curve in 
Figure 6c). This agrees with the findings of Yeates, Constable, and Martens 
(2010), who found that the flux-rope eruption rate increased by a factor of 
eight between 1996 and 1999, comparable to the relative increase of the bipole 
emergence rate, or of the total magnetic energy. 

Note that the rate in 1996 may be underestimated - both in Figure 6 and 
in the earlier simulation - due to a lack of helicity stored in the initial con- 
figuration. Yeates and Mackay (2012) show how the helicity on high-latitude 
PILs depends not only on in situ generation by differential rotation, but also on 
the gradual poleward transport of helicity from lower latitudes. Thus there is a 
time delay for the build up of non-potential structure in the corona. Indeed the 
subsequent Cycle 23 Minimum (2007-2009) shows a higher floor in the eruption 
rate of roughly 0.4 per day. So if the simulation were initialised earlier, we might 
expect a higher eruption rate in 1996. This highlights the importance of long- 
term memory in the coronal magnetic topology. The gradual poleward transport 
of helicity also explains the relatively infrequent eruptions at higher latitudes. 
While such eruptions are present, for example around the "detachment" of the 
polar crowns in 2000, these structures tend to encircle much of the Sun. Once 
helicity is removed in an eruption, it takes time to build again. 



5. Outlook for Space Weather Prediction 

Flux-rope eruptions in the NP model can give us insight into the origin of CMEs, 
although further development is needed before the model can enter the realm of 
operational space-weather prediction. In the first place, the peak eruption rate 
of 1.5 per day in the simulation is substantially lower than the rate of observed 
CMEs, by approximately a factor three when compared with the manually com- 
piled Coordinated Data Analysis Workshops catalogue (CDAW, Yashiro et ai, 
2004), which is based on data from Large Angle and Spectromctric Coronagraph 
(LASCO, Brueckncr et at, 1995). The CDAW rate, divided by three, is shown 
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by the red curve in Figure 6(c). We have filtered out CMEs with apparent width 
less than 15° or greater than 270°, and determined error bars by taking into 
account data gaps (St. Cyr et al., 2000; Ycates, Constable, and Martens, 2010). 

The parameter study of Yeates and Mackay (2009) shows that the eruption 
rate may be increased either by reducing the coronal diffusivity in the simulation 
(equivalent to 774 here), or by increasing the hemispheric twist imbalance of the 
emerging bipolar regions. However, the latter is constrained by observations of fil- 
ament chirality (Yeates, Mackay, and van Ballegooijen, 2008), and for reasonable 
values of these parameters the eruption rate remains too low. 

Many missed eruptions arc likely in active regions. In the present simulations, 
emerging regions are treated as idealised bipoles, and the effect of this simplifica- 
tion is likely to be most pronounced during the early stages of an active region's 
lifetime. Our model was designed to follow the long-term evolution of magnetic 
hclicity, and cannot reproduce multiple eruptions in rapid succession from a 
single active region. Nor are the idealised bipoles a spatially accurate model for 
more complex <5-spot regions. These limitations prevent us from making either 
spatial or temporal predictions of individual CMEs (Yeates et al., 2010b). 

Nevertheless, more detailed modelling of active-region structure and emer- 
gence could in principle be incorporated into the NP model in future. The 
difficulty is that one needs a description of the three-dimensional structure of 
emerging regions. The solution may be either to impose time-dependent electric 
fields on the solar surface (Fan and Gibson, 2007; Fisher et al., 2010), or to 
locally drive the simulation from higher-cadence magnetograms in place of the 
flux-transport model (Mackay, Green, and van Ballegooijen, 2011; Cheung and 
DeRosa, 2012), although these techniques are still under development. 

Wc believe that the time-dependent NP model is worth pursuing because the 
inherent magnetic memory in the corona implies a degree of predictability. While 
flares or CMEs from newly-emerged active regions are unlikely to be predicted 
more than a day or two before the active region appears on the surface, CMEs 
originating from decaying active regions offer scope for longer-term prediction 
using models such as this. Finally, it should be noted that the NP model includes 
only the magnetic field and not other plasma properties. Hence it can predict 
only the initiation of flux-rope ejections, not their subsequent dynamical evolu- 
tion. For individual events, this evolution can be followed in MHD models (e.g., 
Manchester et al, 2004). 



Appendix 

A. Computational Grid 

Our computational grid is divided into latitudinal sub-blocks, each of which has 
uniform spacing in the stretched variables 

x = <f>/A, j/ = -log(tan(0/2))/A, z = log(r/R )/A (16) 

(van Ballegooijen, Priest, and Mackay, 2000), where A is the equatorial grid 
spacing in longitude [</>]. The horizontal cell-sizes arc da; = dy = 1 for the 
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Figure 7. Example of the variable grid with 48 cells at the equator (compared to 192 in the 
actual simulations) and seven sub-blocks (compared to nine). 



equatorial sub-block, and double in each sub-block towards the poles (Figure 7). 
The vertical cell size is dz = 1 for all sub-blocks. This introduction of sub-blocks 
with different spacing counters the problem of grid convergence toward the poles, 
since the horizontal cell area is A 2 r 2 sin 2 Odxdy. The sub-block boundaries in 
latitude are determined so that each grid cell is as large in horizontal area as 
possible, while never exceeding the equatorial cell area A 2 r 2 (Figure 8). With 
a longitudinal resolution of 192 cells at the Equator and 12 at the poles, there 
are nine sub-blocks. (The number 12 is chosen due to the parallel architecture 
used.) The total number of grid cells in (x,y) is 18 936, as compared to 63 744 
for a uniform, single-block grid with unit spacing. 

A complication arising from the variable grid is that the different sub-blocks 
need to communicate ghost values of B$ r and B^ with one another at each 
timestcp. This is analogous to the inter-level communications in Adaptive Mesh 
Refinement (AMR) codes, except that our grid is fixed in time. Restriction 
(from fine to coarse sub-blocks) is the simpler process, for which we use an 
area-weighted average over fine grid cells (Balsara, 2001). Prolongation (from 
coarse to fine sub-blocks) is trickier since the coarse-block solution must be 
interpolated to get the fine-block ghost values at intermediate locations. We 
follow the Taylor expansion method of van der Hoist and Keppens (2007), 
using monotonic van Leer slope estimates (Evans and Hawley, 1988). These 
slopes are also computed throughout the grid and used for slope-limiting in the 
advection terms, in order to prevent spurious oscillations near sharp gradients. 
Our numerical tests indicate that numerical diffusion due to this "upwinding" 
is negligible compared to the physical diffusion D. Finally, after computing the 
update dAo/dt within each sub-block, we replace boundary values on coarser 
sub-blocks with those derived from finer sub-blocks. This is analogous to the 
"flux correction" of Bcrger and Colella (1989). 
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Figure 8. Horizontal cell area (relative to that at the Equator) as a function of latitude, 
where symbols (joined by solid lines) show the variable grid and dashed lines a uniform grid 
with dx = dy = 1 everywhere. 



B. Global Boundary Conditions 

The staggered grid requires ghost-cell values of two components of Bo outside 
each boundary. In longitude the domain is simply periodic. On the photosphere 
r = Rq, we fix ghost cell values of Bog, Bq^ by requiring that vo r = on r = Rq 
in the magncto-frictional model. On the outer boundary r — 2.5R©, we impose 
a radial outflow velocity (see Yeates et at, 2010a) that models the effect of the 
solar wind radially extending field lines, while still allowing horizontal fields to 
escape during flux-rope ejections. The resulting ghost-cell values of B e, B^ do 
not play an important role and are simply set by zero-gradient conditions. On 
the latitudinal boundaries (located at approximately ±89.33° latitude) we need 
ghost cell values for Bg r and Bq^. The latter are simply set to zero, while the 
ghost values of B$ r are chosen to satisfy Stokes' Theorem given the integral of 
Ao^ around the latitudinal boundary. 
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